#!/usr/bin/env python
# -*- coding: utf-8 -*-

"""输出给定栅格的统计值
包括 极大值, 极小值, 均值, 中位数, 标准差 及 直方图
"""

import sys
import cv2
import matplotlib.pyplot as plt
import numpy as np
from osgeo import gdal

def main(input_file):
    gtif = gdal.Open(input_file)

    if gtif is None:
        print(f"Unable to open {input_file}")
        sys.exit(1)

    try:
        srcband = gtif.GetRasterBand(1)
    except RuntimeError as e:
        print(f"No band found")
        print(e)
        sys.exit(1)

    minimum = srcband.GetMinimum()
    maximum = srcband.GetMaximum()
    print(f"Minimum={minimum}, Maximum={maximum}")

    histogram = srcband.GetHistogram(
        min=(minimum-0.5),
        max=(maximum+0.5),
        buckets=int(maximum-minimum)
    )
    # print(f"Histogram={histogram}")
    hist_data = np.array(histogram)
    n, bins, patches = plt.hist(hist_data, normed=True, bins=64)
    plt.xlabel('count')
    plt.ylabel('persent')
    plt.title('DEM Histogram')
    # plt.axis([0, 60, 0, 0.3])
    plt.grid(True)
    plt.savefig('out.png')

    stats = srcband.GetStatistics( True, True )
    print(f"[ STATS ] =  Minimum={stats[0]}, Maximum={stats[1]}, \
Mean={stats[2]}, StdDev={stats[3]}")


if __name__ == '__main__':
    if len(sys.argv) < 2:
        sys.exit(1)

    main(sys.argv[1])
